Draft version January 27, 2012 

Preprint typeset using LAT^X style cmulatoapj v. 5/2/11 



MOLECULAR TRACERS OF TURBULENT SHOCKS IN GIANT MOLECULAR CLOUDS 

A. Pon 1,2 , D. Johnstone 2j , and M. J. Kaufman ;! ' 4 

Draft version January 27, 2012 

ABSTRACT 

Giant molecular clouds contain supersonic turbulence and simulations of magnetohydrodynamic 
turbulence show that these supersonic motions decay in roughly a crossing time, which is less than 
the estimated lifetimes of molecular clouds. Such a situation requires a significant release of energy. 
We run models of C-type shocks propagating into gas with densities around 10 3 cm -3 at velocities 
of a few km s _1 , appropriate for the ambient conditions inside of a molecular cloud, to determine 
which species and transitions dominate the cooling and radiative energy release associated with shock 
cooling of turbulent molecular clouds. We find that these shocks dissipate their energy primarily 
through CO rotational transitions and by compressing pre-existing magnetic fields. We present model 
spectra for these shocks and by combining these models with estimates for the rate of turbulent energy 
dissipation, we show that shock emission should dominate over emission from unshocked gas for mid 
to high rotational transitions (J > 5) of CO. We also find that the turbulent energy dissipation rate is 
roughly equivalent to the cosmic-ray heating rate and that the ambipolar diffusion heating rate may 
be significant, especially in shocked gas. 

Subject headings: ISM: clouds - ISM: molecules - shock waves - stars: Formation - turbulence 
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1. INTRODUCTION 

Molecular line observations of giant molecular clouds 
(GMCs) yield line widths significantly larger than what 
would be expected from therma l motions alone (e.g., 
lLarsorill98lHSoromon et al.lfT987h . These large, nonther- 
mal line widths are generally interpreted as being due to 
supersonic turbulence, with Mach numbers on the order 
oflO (e.g- JZuckerman fc Evansl[l97l iMcKee fe Ostrikerl 
|2007|) . Zeeman splitting measurements of magnetic field 
strengths in molecular clouds show that these supersonic 
motions are on the order of the Alfven speed, which 
suggests that magnetohydrodynamic (M HD) waves may 
play a significant role in molecular clouds (|Crutcherlll999t 
ICrutcher et aJ1l2010h . 

Supersonic, hydrodynamic t urbulence decays on the 
order of a free-fall time (e.g . , iGoldreich fc Kwanl Il974t 
iFieldl IT978L lElmegreenl 119851 ) and thus, maintaining the 
turbulent support of GMCs for their entire lifetimes, 
estimated to be between 2 an d 30 times longer than 
the free fall timescale ( e.g.. iMouschoyias fc Spitzei 
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a sign ificant problem. B ased on theoretical calcula- 
tions (|Arons fc Maxlll975l ). it was believed that MHD 
turbulence would decay an order of magnitude slower 
than hydrodynamic turbulence, thereby preventing the 
dissipation of turbulent energy in GMCs; however, sim- 
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ulations of MHD turbulence show that MHD turbulence 
also d ecays on the order of a free-fall time at the driving 
scale (iGammie fc Ostrikerl Il996t IMac Low et al.l 119981: 
i Stone et al.l 119981: IMac Low! Il999t iPadoan fc Nordhmdl 
ll999HOstriker et al.ll200"lll 

In MHD turbulence simulations, turbulent energy is 
dissipated via numerical viscosity and artificial viscosity 
in shock fronts. Under the assumption that the dissi- 
pated turbulent energy is lost as heat and rapidly radi- 
ated away, many MHD simulations are run with isother- 
mal equations of state and thus, these simulations do 
not explicitly fo llow where the dissipated turbulent en - 
ergy goes (e.g. , IStone et al.l 119981: ISmith et al.ll2000bl) . 
Bas u fc Muralil (l20o¥ made a first attempt to compare 
the CO J = 1 — » luminosities of molecular clouds to 
what they predicted would be seen from molecular clouds 
based upon simple energetic arguments. Since then, how- 
ever, little progress has been made in determining where 
this turbulent energy goes and whether there are any 
observational signatures of this dissipated energy. 

Shocks increase the temperature and density of the 
shocked gas, which, in turn, can substantially alter the 
chemistry of the gas and the emissio n coming from the 
gas (e.g.. Kaufman fc Neufeldlll996bl ia|). thereby poten- 
tially providing a distinct signature and tracer of tur- 
bulent energy dissipation via shocks. For typical turbu- 
lent velocities and magnetic field strengths of molecular 
clouds, the magnetic field is capable of transmitting in- 
formation about the presence of a shock to ions upstream 
of the shock front. This eliminates discontinuities in gas 
properties across the shock front and spreads out the 
thickness of the shock. In turn, this leads to lower tem- 
peratures in the shocked gas and prevents molecules from 
being dissociated. Such a shock is referred to as a con- 
ti nuous, or C-typ e , shock and is described in more detail 
in iMullanl (fl97l . iDralnel ([1980), and iDraine fc McKed 
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iKaufman fc Neufeldl (|1996aD . propagating into molecu- 
lar gas with densities around 10 3 cm -3 at velocities of 
a few km s _1 , appropriate for the ambient conditions 
inside of a molecular cloud, to determine which species 
and transitions dominate the cooling and radiative en- 
ergy release associated with shock cooling in turbulent 
molecular clouds. The shock velocities modeled are on 
the order of the typical turbulent velocity of molecular 
clouds, which are much lower than the velocities of pro- 
tostellar outfl ows that have been the target of previous 
studies (e.g.. iChernoff et aL Il982t iTimmermannl 119961 : 
IKaufman fc Neufeldl Il996bl ]al). These shock models are 
combined with estimates for the rate of turbulent energy 
dissipation in molecular clouds to predict the integrated 
intensities of various shock excited lines coming from an 
entire molecular cloud and these integrated intensities 
are then compared to those fr om photodissociation r e- 
gion (PDR) models based upon Kaufma n et al.1 £1999). 

Typical scaling relations of GMCs are presented in 
Section 12.11 and the turbulent energy dissipation rate of 
molecular clouds is derived in Section l2~2l The shock and 
PDR models used in this paper are described in Sections 
12. 3112. 51 and the results of these models are presented in 
Section [3] In Section |H the implications of these results 
are discussed, and in Section [SJ the rate of turbulent dis- 
sipation is compared to other known heating mechanisms 
in molecular clouds. Finally, our findings are summarized 
in Section O 

2. SETUP 

2.1. Scaling Relations of GMCs 

Correlations between the size, density, and line-of-sight 
velocity dispersion of GMCs, as determined through CO 
observations, are well known and are coll e ctively referred 
to as Larson's laws (e.g., ILarsonl Il981t ISolomon et al.l 



119871: IHever fc Brvmtl 120041) . Th e best fitting scaling re 
lations found by ISolomon et al.1 (|1987l ) are: 



cr = 0.72(i?/pc) a5 km s" 1 , 



134( J R/pc)" 1 M Q pc 
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(1) 

(2) 



where R is the effective radius (the radius of a spher- 
ical cloud with the same projected surface area as the 
observed cloud), a is the one-dimensional velocity dis- 
persion (which we assume is equal to the observed line- 
of-sight velocity dispersion), and p is the average density. 
These relations suggest that a cloud with a mean density 
of 10 3 cm" 3 has a radius of approximately 2 pc, a mass 
of 2000 M , a total molecular hydrogen column density 
of 1.2 x 10 22 cm" 2 (corresponding to a visual extinction 
of 12 through the entire cloud), and a one-dimensional 
velocity dispersion of about 1 km s . 

The size- velocity relationship is fairly well established, 
although there is some evidence that the velocity dis- 
persion of a molecular cloud may also depend upon the 
column density of that cloud (|Hever et al.ll2009h ■ The va- 
lidity of the size-density relation, however, is much less 
certain, as the observed relationship may be only due 
to the limited dynamical range of cur rent observations 
(|Ballesteros-Paredes fc Mac Lowll2002D . For the shock 
models used in this paper, Larson's laws are only used to 
confirm that the simulated parameter range roughly cor- 
responds to the properties of observed molecular clouds. 



Neither part of Larson's laws is used to calculate the in- 
t egrated intensity of th e shock emission. 

ISolomon et all (|1987f) also found a correlation between 
the velocity dispersion and 12 CO J = 1 —> total lumi- 
nosity , Li_»p, of mole c ular c louds. In units of K km s" 1 
pc" 2 , ISolomon et~aT1 (|1987ft found that the CO 1 ->• 
luminosity is: 

L^ = 130a 5 . (3) 

The relationship between the magnetic field strength 
in a molecular cloud, B, and the number density of hy- 
drogen nuclei, uh , is often expressed in the form 



B = bn k H pG, 



(4) 



where b and k are the fitting parameters. A value of k 
= 0. 5 corresponds to a cons tant magnetic energy den- 
sity (|McKee fc Ostrik"erll2007l ) a nd is expected from am- 
bipol ar diffusion collapse models dFiedler fc Mouschoviasl 
1993). A k = 0.5 relation is also expected if the tur- 
bulent velo city in a cloud is always roughly the Alfven 
speed (e.g., ICrutcherl[T9"9 9'). A value of k = 2/3, how- 
ever, is predicted if magnetic fields are unimportant and 
a molecular cloud is able t o maintain a roughly spherical 
shape during its col lapse (jMestell 119661 : iCrutcherl 119991 : 
ICrutcher et ffllMToh . 

ICrutcherl (|1999l ) compiled Zeeman splitting observa- 
tions and found b = 0.95 and k = 0.5 (Mc Kee fc Os trikcr 
120071) The MHD simulations of IPadoan fc Nordlundl 
( 1999) exhibit a k — 0.4 relation and the relation b = 1 , 
k = 0.5 is commonly adop ted (e.g., IDraine "eTaIlll98l 
IKaufman fc Neufeldl 1 1996bT iaT). Recently. ICrutcher et al.1 
(|2010f ) have examined all of the available Zee man split- 
ting o bservations, including those used by ICrutcherl 
(1999), and found that the best fit for the maximum ob- 
served line-of-sight magnetic field strength comes from 
the relation: 



Br, 




n H < 300 cm" 3 



n H > 300 cm 



(5) 



For densities greater than 300 cm" 3 , the above relation 
corres ponds to k = 2/3 and b = 0.25. ICrutcher et al.1 
(2010) also note that their data arc consistent with hav- 
ing line-of-sight magnetic field strengths down to essen- 
tially zero. 

The maximum magne tic field strengths that were fit 
by ICrutcher et all (|2010f ) most likely correspond to cases 
where the magnetic field is highly aligned with the line- 
of-sight, such that the full magnetic field strength is mea- 
sured. The average magnetic field strength along any 
random direction will thus be only half of the strength 
given by the above relatio n. For a cloud with an H2 
density of 10 3 cm" 3 , the ICrutcher et~aTl 1)20100 rela- 
tion therefore predicts that the average magnetic field 
strength along any random direction is 17 pG. Intrinsic 
scatter in the magnetic field strength between different 
clouds will also likely further reduce the average mag- 
netic field strength. 

2.2. Turbulent Energy Dissipation Rate 

The turbulent energy density of a molecular cloud is 
approximately 



Eturb — ^ P 



(6) 
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Following the discussion in B asu fe Muralil (|2001l) , the 
mean turbulent energy dissipation rate per volume can 
be written as T tur b = Eturb/td, where td is the dissipation 
timescale. We define the flow crossing time of the cloud 
as t c — 2R/o~ and introduce the ratio of the dissipation 
time to the flow crossing time as a new parameter: k — 
td/t c . The turbulent dissipation rate per volume is thus 



is 



turb 



3/9 cr 3 
AkR' 



(7) 



As shown in lBasu fc M urali (200f ), rather than writing 
the turbulent energy dissipation rate in terms of k, the 
dissipation rate can be expressed in terms of the driving 
scale of the turbulence, A: 



r 



turb 



pa" 



(8) 



where rj is a dimensionless parameter that is a function of 
the density, velocity dispersion, and driving wavelength. 
Comparing Equation ([7]) to Equation © gives the rela- 
tion 

K= . 9 

4rjR V ' 

Periodic box simulations of MHD turbulence have 
found that for a variety of initial conditions, 77 has 
a value between 0.5 and 4 (IGammie fe Ostriker l Il99l 
IStone et all H998I: IMac Low! 119991: IQstriker et al.l 12001). 
Unfortunately, there is no clear consensus on what scale 
turbulence is driven on. Protostellar outflows, which 
drive turbulence on small scales, appear to have enough 
energy to drive turbulence in active star forming region s 
(|Quillen et al.ll2005t ICurtis et al.l[2010l: lArce et al.ll2010h . 
It is, however, unclear whether outflows are capable of 
driving turbulence a c ross an entire m olecular complex 
(|Baneriee et al.ll2007t lArce et al.lfeOlOl ). Studies of den- 
sity and velocity structure in molecular clouds find that 
the observed structures are only consistent with driv- 
ing at size scales at, or above, the size of the cloud (e.g., 
Ossenkopf fe Mac | Lowl | 2002HBrundl2003HHever fe Bruntl 
20041: iBrunt et al.ll2009l; iPadoan et al.ll2009H . Supersonic 



turbulence has also been observed in the Polaris Fl are, 
which is devoid of any protostars ()Andre et al.ll2010h . 

For the remainder of this paper, a n value of one will be 
adopted, for which the turbulent dissipation timescale is 
equal to the flow crossing time of the cloud. Via Equation 
(0) and the numerical factors for 77, this corresponds to 
a turbulent driving scale on the order of the size of a 
molecular cloud. 

Equation (|7|) is a general result that can be applied to 
any cloud, given that the characteristic radius, density, 
and velocity dispersion are known. For this paper, the 
simplifying assumption that clouds are spherical will be 
made, such that the total turbulent energy dissipation 
rate is 



I'turb — 



J turb " 



3/9 cr 3 4tt i? 3 
AkR 3 
ir pa 3 R 2 



(10) 
(11) 



If all of the dissipated turbulent energy is radiated 
away, the corresponding total integrated intensity, Iturb, 



Iturb 



urb 



4tt 2 R 2 ' 
3 



I (u 



pcr_ 

47T K 



(12) 
(13) 



This integrated intensity is independent of the size of 
the molecular cloud. For this paper, a mean mass per 
particle of 4.6 x 10~ 24 g, or about 2.77 amu, is adopted, 
and with this value, the above equation becomes 



Iturb = 3.66 x 10 



(l0 3 cm- 3 ) (l kms- 1 ) 



10 3 cm 
erg s -1 cm -2 steradian -1 



(14) 



J turb = 8.60 x 1Q- 1S k- 1 ( 



10 3 cm 3 / V 1 km s 



erg s 1 cm 2 arcsec 2 . 

2.3. Shock Code 



(15) 



To determine which species and transitions dominate 
the cooling and radiative energy release associated with 
shock cooling of turbulent molecular clouds, we run 
models of C-type shocks based upon the models of 
iKaufman fe Neufeldl (|1996al ) with initial conditions cor- 
responding to that expected for roughly one parsec sized 
molecular clouds. 

The code used first calculates the temperature, den- 
sity, chemical abundance, and velocity profiles of a 
C-type shock. To do so, it calculates the cool- 
ing rates for rotatio nal and vibrational t ransit ions of 
H 2 0, H 2 , and CO (jNeufeld fc Kaufman! fl993T) : colli- 
sions between the neutral gas and cooler dust grains 
(Hollenbach & McKec 1989) : and H2 dissoc i ative cooling 
(|Lepp fc Shulll 119831 IHollenbach fc McKeel [l989Tl . The 
freezing out of molecules on dust grains is, however, 
not calculated by the code as the freezeout timescale 
is expected to be much longer than the shock cooling 
timescale. 

Once the shock structure is determined, the code then 
calculates the integrated intensities of each molecular 
transition of interest by solving the partial differential 
equations for the line emission at each point and then 
integrating the emission over the entire shock profile. 
This extra step of determining individual line strengths, 
rather than just determining an overall cooling rate for 
a particular molecule, is only used for CO. In this later 
step, the code only includ es rotational line emi s sion fo r 
gas down to 10 K, unlike in IKaufman fc Neufeldl (|1996a|) , 
where emission is only included from gas above 50K. No 
such temperature limitation is used when calculating the 
overall cooling rates for each molecule in the first half 
of the code. For a mo re detailed description of how this 
code works, please see Kaufma n fc Neufeldl (|1996al) . 

By equating the total kinetic energy dissipated by 
shocks to the turbulent energy dissipation rate of a 
molecular cloud, given by Equation (|11|) . we scale our 
shock models to predict the expected integrated intensi- 
ties from each CO rotational line. That is, the integrated 
intensity of each line is set to the appropriate fraction of 
the integrated intensity given by Equation (|T5|) . It is as- 
sumed that the shock emission is coming from a region 
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larger than the size of the beam and each line is scaled 
equally under the assumption that all of the lines are 
optically thin. The effect of the lower lying lines being 
optically thick is discussed further in Section l4~6l 



Table 1 

Shock Model Properties 



2.4. Shock Code Parameters 

For each shock model, the same, roughly solar, chem i- 
cal composition as used in lKaufman fc Neufeldl (|1996bl |aj) 
is used. In particular, the initial CO number density is 
set to be 1.2 x 10~ 4 times that of the H nuclei number 
density and the initial H2O abundance is set to 10 -7 . As 
shown in Section [2.1l a one parsec sized molecular cloud 
is expected to have a density of approximately 10 3 cm 3 . 
Thus, the initial H 2 density is set to be either 10 2 5 , 10 3 , 
or 10 35 cm -3 in these models. 

If the velocity distribution of gas particles in a molec- 
ular cloud is Gaussian in every direction, with a one- 
dimensional velocity dispersion of a, then the distribu- 
tion of relative velocities between two gas particles in 
the cloud will also be Gaussian in every direction with 
a one-dimensional velocity dispersion of y/2cr. Since the 
energy dissipation rate of a shock scales with the third 
power of the shock speed, the mean speed at which en- 
ergy is dissipated is the cube root of the mean cubed ve- 
locity difference between two gas particles, < Av 3 > 1 / 3 , 
which is roughly 2.4c The shock velocity at which the 
peak energy dissipation rate occurs is slightly higher, 
approximately 3.2er. Thus, the characteristic shock ve- 
locity in a molecular cloud with a one-dimensional ve- 
locity dispersion of 1 km s _1 , consistent with the size- 
velocity relation for a radius of 1 pc, is on the order of 
2-3 km s _1 . For the remainder of this paper, we assume 
that the one-dimensional velocity dispersion is a factor 
of 3.2 smaller than the shock velocity. The larger con- 
version factor of the two mentioned above is chosen so 
that the corresponding velocity dispersions, and thus the 
shock-integrated intensities calculated in Section [3~Tl are 
smaller. 

Models with shock velocities of 2 and 3 km s" 1 are 
computed. For a temperature of 10 K, these velocities 
correspond to Mach numbers of 12 and 17, respectively. 
While these velocities are appropriate for turbulent mo- 
tions in a molecular cloud, they are much lower than 
the typical velocities of protostellar outflows and winds. 
Such higher velocity flows have been modeled extensively 
in the past and give rise to significantly higher p ost shock 
temperatures (e.g.. Kaufman fc Ncufcld 1996b|a]). 

The strength of the magnetic field parallel to the shock 
front is initialized using the parameterizations k = 0.5 
and b = 0.1 or 0.3, where b and k are as defined in 
Equation (j4]). The component of the magnetic field per- 
pendicular to the shock front is always set to zero, as 
this component has no effect on the shock structure in 
our steady state, plane parallel models. Thus, the initial 
magnetic field strength ranges from 3 /iG to 24 fiG in 
the different models. For a weaker field parallel to the 
shock front, the shock thickness is smaller and the en- 
ergy released in line radiation is relatively larger. This 
is why magnetic field strengths that are slightly lower 
than, although still generally consistent with, the average 
line-of-sig ht magnetic field stren gth given by the scaling 
relation of iCrutcher et all (|2010[) have been chosen. 
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Note. - 


Column 


1 shows the 


model 


name, 


while Columns 



2 and 3 show the logarithm of the initial density and shock 
velocity of each model respectively. Column 4 represents the 
magnetic b parameter, as defined in Equation |4j, and Column 
5 shows the resulting initial magnetic field strength. Columns 6 
and 7 show the Mach number and Alfvenic Mach number of the 
models, respectively. 



The Alfvcn speed is 



D 



VA 



(16) 



For our shock models, the Alfvenic Mach number, given 
by M A = v shock /v A , ranges from 4 to 16. 

Twelve shock models are run, one for each combination 
of initial density, shock velocity, and magnetic field b 
parameter. Table Q] gives the shock velocity, magnetic 
field strength, Mach number, and Alfvenic Mach number 
for each model. A naming convention of nWXvYbZ is 
adopted, where W.X is the logarithm of the initial H2 
number density in cm~ 3 , Y is the shock velocity in km 
s , and Z is the magnetic b parameter. 

2.5. PDR Model 

The shocked gas in a molecular cloud is not the only 
source of molecular line emission. The cool, well-shielded 
gas and the warmer gas in the PDR at the cloud's 
surface, which is exposed to the interstellar radiation 
field (ISRF), will also contribute emission. To model 
thi s emission from unsho cked gas, PDR models based 
on Kauf man et all (119991) are used. As suggested by 
iKaufman et al.l ( 19991) . these plane parallel models are 
adjusted for spherical geometry by using the equation 



L = I A7rj(N)r 2 dr, 



(17) 



where r is the radial distance from the center of the cloud 
and j (N) is th e emissivity at a c olumn N from the surface 
of the cloud. [Kaufma n et al.l ()1999[ ) estimate that this 
procedure produces results that are within a factor of 1.5 
from intrinsically spherical PDR models. Furthermore, 
for any optically thin line, the resulting integrated inten- 
sity is doubled to account for photons originally emitted 
radially inward. 

The PDR models used have ISRFs of 3 Habing, where 
the average far ultraviolet ISRF i n free space is 1.7 
Habing or 2.7 x 10" 3 erg cm" 2 s" 1 (jTielensll2005h . and 
microturbulent Doppler line widths of 1.5 km s , sim- 
ilar to the velocity dispersions of the shock models. It 
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is assumed that the PDR emission fills the beam. These 
PDR models do not take into account the freezing out of 
CO onto dust grains. 

A density of 10 3 H n udei cm~ 3 is used for all of the 
comparison PDR models, which is comparable to the me- 
dian initial density in the shock models. We believe that 
this is an appropriate comparison density for the 10 3 ' 5 
cm~ 3 shock models because a density gradient should 
be present within realistic molecular clouds, with the 
density decreasing toward the periphery of the cloud, in 
order for the clouds to remain in pressure equilibrium. 
Thus, it is expected that the warm outer layers of molec- 
ular clouds, from which most of the PDR emission comes 
from, are at lower densities than the bulk of the cool, CO- 
rich gas in the interior of the clouds from which most of 
the shock emission originates. In PDR models with den- 
sities below 10 3 Hnuciei cm -3 , CO only forms in the gas 
phase once the gas has cooled significantly, such that 
there is almost no emission in the mid to high J lines 
of CO. To be more conservative in our findings of when 
shock emission is stronger than PDR emission, we pre- 
fer using PDR models with densities of 10 3 Hnuciei cm~ 3 
for comparison with the 10 3 and 10 25 cm~ 3 shock mod- 
els, even though these PDR models may over predict the 
CO emission from 10 25 cm -3 gas. The results should be 
roughly consistent for the 10 3 cm -3 shock models. 

For the comparison PDR models, the size of the molec- 
ular cloud must also be known in order to know at 
what Ay to cut the model. For each shock model, the 
ISolomon et al.l (|1987T ) size- velocity relation is used to de- 
termine an appropriate, typical size for a cloud and then 
the CO column density of that cloud is determined un- 
der the assumptions that the cloud i s spherical and has 
a typical CO abundance of 2 x 1(T 4 ([Glover et al.ll2010T> 
throughout the entire cloud. This CO column density is 
then used to determine the appropriate depth of the com- 
parison PDR model such that the PDR model has the 
same CO column density. The depths of the PDR models 
are chosen based upon a CO column density, rather than 
upon a hydrogen nuclei column density, because the ini- 
tial chemical abundances used for the shock models are 
consistent with gas in which CO has already formed in 
the gas phase and because CO cooling dominates the en- 
ergy b udge t of the shock models, as described later in 
Section 13.11 

Two shock models, n35v3bl and n35v3b3, require CO 
column densities larger than the total CO column present 
at the maximum depth of the 10 3 cm~ 3 PDR model. The 
contribution from the most deeply embedded layers of 
this PDR model to the total emergent flux does, however, 
drop to negligible values for all of the CO lines. This is 
because the lower lines are optically thick and the gas 
temperature at high Ay is too low for any significant 
emission in the higher lines. Thus, we use the full extent 
of the 10 3 cm~ 3 PDR model as the comparison model 
for these two shock models and we do not believe that 
this failure to exactly match the CO column densities of 
these two models significantly affects our results. 

2.6. Empirical CO 1^0 Luminosities 

As described in Section [21] ISolomon et al.l ([19871 ) 
found an empirical correlation between the velocity dis- 
persion and 12 CO J = l->0 total luminosity of molecular 
clouds. For the velocity dispersion corresponding to each 



of the shock mo dels, the 1 — » integr ated intensity ex- 
pected from this ISolomon et al.l (119871) relation is calcu- 
lated using a cloud radius from the ISolomon 
size-velocity relation and the assumption that molecular 
clouds are spherical. While we have used the size- velocity 
relationship in determining the comparison PDR spectra 
and in calculating the empirically expected 1 — >■ inte- 
grated intensity, we re-emphasize that this relationship 
is not used at any time in the calculation of the expected 
shock spectra. 

3. RESULTS 
3.1. Shock Line Emission 

In the upper panels of Figure [1] the neutral velocity, 
ion velocity, and density profiles of models n30v2bl and 
n30v3bl are shown. As typical in MHD, C-type shocks, 
the density and velocity profiles show no sharp discon- 
tinuities and the ion velocity decreases before the neu- 
tral velocity does. Due to mass conservation, the den- 
sity and neutral velocity are inversely related to each 
other, and thus, the maximum density reached in model 
n30v3bl is larger than that reached in model n30v2bl. 
The magnetic field strength and ion velocity are similarly 
inversely correlated. 

The lower panels of Figure[T]show the temperature pro- 
files of models n30v2bl and n30v3bl as well as the cool- 
ing profiles due to CO, H2, and H2O lines and gas-grain 
coupling. During the initial stages of a shock, where the 
gas temperature is still increasing, the rate of CO cooling 
increases in close tandem to the increase in temperature. 
After the gas has reached its peak temperature, the CO 
cooling rate remains higher than it was at the same tem- 
perature earlier in the shock. This is due to the higher 
gas densities in these later regions of the shock that al- 
low for more efficient population of higher J CO states 
and thus, more effective CO cooling. The CO cooling 
rate is the least temperature sensitive of any of the plot- 
ted cooling terms, whereas the H2 cooling rate shows 
a strong dependence on temperature. The H2 cooling 
rate is strongly peaked around the temperature peak and 
shows the most significant change between the two shock 
models. While the gas-grain cooling rate is temperature 
dependent, it is also clearly larger at higher densities as 
the gas-grain cooling curves are skewed toward the higher 
density sides of the shocks. High frequency noise, which 
is likely numerical in nature, appears toward the end of 
both models in the CO cooling rate and thus, a boxcar 
smoothing algorithm has been applied to the CO cooling 
rates for distances larger than 0.06 pc in model n30v2bl 
and for distances larger than 0.04 pc in model n30v3bl. 
This noise does not significantly affect our results as it 
only occurs over very limited spatial scales and occurs 
only when the cooling rates have decreased significantly, 
such that the noise does not significantly affect the total 
cooling rate. Furthermore, this noise only occurs when 
the temperature has dropped below 10 K, at which point 
the lack of cosmic-ray heating in our models becomes im- 
portant (see Section |4~T|) . The other shock models are 
qualitatively similar to the ones shown in Figure Q] and 
thus, are not shown. The temperatures, densities, and 
timescales of the shocks are too low for any significant 
chemical changes to occur within the gas in any of the 
models and thus, the small changes in chemical abun- 
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dance across the shock models are also not shown. 
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Figure 1. Various profiles of models n30v2bl and n30v3bl. The 
top row shows density, neutral velocity, and ion velocity profiles as 
the solid (black), dotted (blue), and dashed (red) lines, respectively. 
The velocity axis is given on the left-hand border while the density 
axis is given on the right-hand border. The bottom row shows 
temperature profiles as the solid (black) lines and cooling profiles 
due to CO, H2, gas-grain interactions, and H2O as the dotted 
(blue), dashed (red), dash-dotted (green), and dash-triple-dotted 
(yellow) lines, respectively. The cooling rate axis is given on the left 
border and the temperature axis is given on the right border. The 
CO cooling profiles have been boxcar smoothed beyond a distance 
of 0.06 pc in model n30v2bl and past 0.04 pc in model n30v3bl 
due to the presence of high frequency noise. This noise is likely 
numerical in nature an d sho uld not significantly affect our results 
(see the text in Section l3.ll . The left-hand column shows profiles 
of the n30v2bl model and the right-hand column shows profiles 
of the n30v3bl model. The x-axes of all four boxes are the same 
and the y-axis scaling is the same for both models. Please see the 
online version for a color version of this figure. 

The dominant molecular coolant in all of these slow 
shock models is 12 CO, with 40%-80% of the dissipated 
energy going into 12 CO rotational lines. A significant 
fraction, 15%-60%, of the dissipated energy is not ra- 
diated away, but rather, goes toward compressing the 
magnetic field. In the models with the weakest shocks, 
those with b = 0.3 and a shock velocity of 2 km s _1 , 
the conversion of kinetic energy into magnetic energy is 
the most significant mechanism for dissipating kinetic en- 
ergy. Molecular hydrogen rotational lines are the second 
most effective molecular coolant, but dissipate less than 
1% of the shock energy in all but the models with b = 
0.1 and a shock velocity of 3 km s" 1 , which are the mod- 
els with the strongest shocks. In these stronger shock 
models, H2 lines account for between 7% and 21% of the 
energy dissipated, with H2 cooling being more important 
at lower densities. All other cooling mechanisms are very 
minor in these shock models. A summary of where the 
energy goes in each model is given in Table [5J It should 
be noted that the sums of the cooling functions do not 
exactly equal the total kinetic energy dissipation rates. 



Table 2 

Sources of Energy Dissipation in the Shock 
Models 



Model 
(1) 


Enn 
(2) 


Eb 
(3) 


Eh 
(4) 


Edust 

(5) 


(ef 


IlinJ V £i U _1_ 


76 


23 


8 


2 


0.02 


n25v3bl 


61 


17 


21 


2 


0.02 


n25v2b3 


12 


56 


<0.1 


<0.1 


0.01 


n25v3b3 


54 


4-1 


0.9 


1 


0.01 


n30v2bl 


76 


26 


0.2 


2 


0.02 


n30v3bl 


68 


18 


11 


3 


0.03 


n30v2b3 


43 


54 


<0.1 


0.2 


0.01 


n30v3b3 


55 


41 


0.3 


2 


0.02 


n35v2bl 


74 


25 


0.1 


3 


0.04 


n35v3bl 


72 


19 


7 


4 


0.04 


n35v2b3 


41 


53 


<0.1 


1 


0.03 


n35v3b3 


53 


42 


0.1 


3 


0.03 


Note. 


— Column 1 


shows the 


model 


names. 



Columns 2-6, respectively, list the percentage of the 
kinetic energy of the shock that is dissipated via 
CO rotational lines, increasing the magnetic field 
strength, H2 lines, gas-grain collisions, and H2O 
lines. 

The total cooling rate, however, is never more than 5% 
discrepant from the kinetic energy dissipation rate. We 
believe that this discrepancy arises from difficulties in ex- 
tending our shock cooling functions to low temperature 
but do not believe that this small discrepancy signifi- 
cantly affects our results. 

The 12 modeled shocks are relatively weak shocks 
and produce density enhancements of at most a fac- 
tor of ^20, of the order of the Mach number. Such 
compressions are still smaller than the density contrast 
between the ambient material in molecular clouds and 
dense cores, which have densit ies of 10 5 cm" 3 or above 
(e.g.. iDi Francesco et alj|2007f) . The maximum temper- 
ature in the shocked gas varies significantly, with the 
stronger shock models achieving maximum temperatures 
of approximately 150 K and the weaker shock models not 
even warming up to 20 K. The maximum density and 
temperature reached in each model is given in Table |3l 

The cooling length of each shock is taken to be the 
full width at quarter maximum of the total cooling func- 
tion profile. The cooling lengths range from 0.01 pc to 
0.35 pc, with the high magnetic field strength and low- 
density models having the largest cooling lengths. The 
corresponding cooling timescales range from 5 x 10 3 years 
to 3 x 10 5 years, with the longer cooling timescales cor- 
responding to larger cooling lengths. 

The volume filling factor of shocked gas in a molecular 
cloud, //, can be calculated from 

» » r turb d CO ol f-\ a\ 

ff= AE k ' (18) 

where T tur f, is the turbulent energy dissipation rate per 
volume (given by Equation (J7J), d coo i is the cooling 
length of the shock, and AE k is the kinetic energy dis- 
sipated per shock front area. Both d coo i and AEk are 
calculated by the shock code but the cloud radius is re- 
quired to calculate Tturb- For this filling factor calcula- 
tion, we use the relatively well-established size-velocity 
relation of Larson's laws, Equation (JT|), to determine the 
appropriate cloud radius for each shock model. We reit- 
erate, however, that the integrated intensities presented 
in Figures [2]|4] and in Table [4] are derived independently 
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Table 3 

Shock Model Properties 



-17- Model n25v2b1 



Model 


log(n maa; ) 


'-L j-ficix 


tcool 


dcool 


ff 


Rb 


Rad 




(cm -3 ) 


(K) 


(10^ years) 


(pc) 


(%) 






(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 

V ) 


(8) 


niovzbl 


3.7 


56 


6.4 


0.06 


0.38 


14 


54 


IlZitJ VOU1 


3.8 


1 15 


3.0 


0.04 


0.11 


22 


104 


n25v2b3 


3.1 


11 


28.5 


0.35 


2.15 


4 


10 


n25v3b3 


3.3 


60 


9.8 


0.17 


0.46 


7 


19 


n30v2bl 


4.1 


54 


2.6 


0.03 


0.16 


11 


57 


n30v3bl 


4.3 


151 


1.3 


0.02 


0.05 


22 


106 


n30v2b3 


3.6 


13 


11.8 


0.14 


0.87 


4 


9 


n30v3b3 


3.S 


60 


3.9 


0.07 


0.18 


7 


18 


n35v2bl 


4.6 


53 


1.1 


0.01 


0.07 


14 


55 


n35v3bl 


4.7 


157 


0.5 


0.01 


0.02 


21 


108 


n35v2b3 


4.1 


17 


1.9 


0.06 


0.36 


4 


9 


n35v3b3 


4.3 


61 


1.7 


0.03 


0.08 


7 


18 


Note. - 


Column 1 shows the model names 


while Columns 2 


and 3 show 



the logarithm of the maximum density reached and the maximum temperature 
reached, respectively. Column 4 shows the cooling time of the shocked gas 
and Column 5 shows the cooling length. The volume filling factor of shocked 
gas for a cloud compatible with the size-velocity relationship of molecular 
clouds is given in Column 6. Columns 7 and 8 show the factors by which 
the magnetic field strength and the ambipolar diffusion volume heating rate 
increase between the initial and shocked gas. 
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from any part of Larson's laws. The volume filling factor 
of shocked gas is always between 0.02% and 0.5% of the 
cloud volume, except for the three weakest shock models, 
where the volume filling factor becomes as large as 2%. 
The cooling times, cooling lengths, and filling factors for 
all 12 models are given in Table [3J 

Figure [5] shows the integrated intensities of CO rota- 
tional transitions as calculated from the shock models 
with densities of 10 25 cm~ 3 . The CO spectra from the 
corresponding comparison PDR models, as well as the es- 
timates for the J = 1 — > integrated intensities from the 
iSolomon et al.l (|1987l ) scaling relation, are also shown in 
Figure [5J Figures [3] and @] show the shock spectra from 
the models with densities of 10 3 cm -3 and 10 3 5 cm~ 3 
resp ectively, as wel l as th e corresponding PDR spectra 
and ISolomon et al.1 ()1987f ) J — 1 — >• integrated intensi- 
ties. All integrated intensities are given in units of erg 
s _1 cm -2 arcsec -2 . The y-axis ranges in Figures [21 G2 
and S] are identical to facilitate comparisons between the 
different shock models. As explained in further detail in 
Section |4~T1 CO spectra have not been calculated for any 
shock models with b = 0.3 and a shock velocity of 2 km 
b- 1 . 

In all of the cases considered, the PDR emission is sig- 
nificantly stronger, by approximately an order of mag- 
nitude, than the shock emi ssion in the thr e e lowe st CO 
rotational transitions. The ISolomon et al.l (|1987l) 12 CO 
J = 1 — >• integrated intensity is also larger than the 
predicted shock J = 1 — > integrated intensity in all 
of the models. On the other hand, in all of the mod- 
els, the shock emission is stronger for all of the high J 
transitions. In the strongest, high-density shock models, 
the shock-integrated intensity becomes larger than the 
PDR-intcgrated intensity in the J = 5 — > 4 line. The 
PDR emission in the higher J lines drops much more 
rapidly than the shock emission and the shock emission 
is an order of magnitude stronger in the J — 7 — > 6 line 
in most models. 

While the shock code does not calculate the individual 
strengths of the different H2 rotational transitions, the 



Figure 2. Integrated intensities of various 1 CO rotational tran- 
sitions for shock models with densities of 10 2 ' 5 cm~ 3 in units of 
erg s — 1 cm -2 arcsec -2 . The shock velocity and magnetic field b 
parameter used for each shock model are given on the top and left 
of the grid, respectively, while the model name is given in the top 
left hand corner of each box. The green (lightest) lines show the 
ISolomon etaLl fl987T ) 12 CO J = 1 -> line strengths, the blue 
(darkest) lines show the CO shock spectra, and the red (medium) 
lines show the CO PDR spectra. The ten lowest rotational tran- 
sitions of CO are labeled in the lower right grid panel. Note how 
the shock spectra dominate over the PDR spectra for high J tran- 
sitions. 

lowest rotational levels should be in local thermodynamic 
equilibrium (LTE) at the modeled densities. Thus, it 
is assumed that all of the H2 line ratios are given by 
their LTE values. At 50 K, the ratio of S(l) to S(0) is 
only 3.3 x 10 -3 , but at 150 K, this ratio increases to 2.3. 
The S(2) and higher lines are negligible in comparison to 
the S(l) and S(0) lines at these temperatures. For the 
three models with significant H2 emission, the expected 
integrated intensities of the S(l) and S(0) lines are given 
in Table |H 

While the PDR models used do not include the 
expected line strengt hs for H2 rotational transitions, 
Kaufman etHI (f2006l ) found that the S(l) and S(0) H 2 
lines have integrated intensities of 7x 10 -20 erg s" 1 cm~ 2 
arcsecond -2 and 2 x 10 -19 erg s _1 cm -2 arcsecond -2 , re- 
spectively, in a PDR with a density of 10 3 cm -3 and an 
ISRF of 3 Habing. Therefore, the PDR H 2 -integrated 
intensities are expected to be comparable to or slightly 
lower than the shock H2-integrated intensities. 

4. DISCUSSION 

4.1. CO Lines 

As mentioned in Section l3T| no CO spectra are calcu- 
lated for the models with b = 0.3 and a shock velocity 
of 2 km s . This is because the maximum post shock 
temperatures in these three models, 11, 13, and 17 K, are 
on the order of the te mperature expected from cosmic- 
ray h eating alone (e.g., iGoldsmithl 120011 : iPan fc Padoanl 
2009). The shock code used does not contain any addi- 
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Figure 3. Integrated intensities of various 1 CO rotational tran- 
sitions for shock models with densities of 10 3 cm -3 in units of erg 
s _1 cm -2 arcsec -2 . The shock velocity and magnetic field b pa- 
rameter used for each shock model are given on the top and left 
of the grid, respectively, while the model name is given in the top 
left hand corner of each box. The green (lightest) lines show the 
ISolomon et~aH I|1987I 1 12 CO J = 1 -> line strengths, the blue 
(darkest) lines show the CO shock spectra, and the red (medium) 
lines show the CO PDR spectra. The ten lowest rotational tran- 
sitions of CO are labeled in the lower right grid panel. Note how 
the shock spectra dominate over the PDR spectra for high J tran- 
sitions. 



Figure 4. Integrated intensities of various 1 CO rotational tran- 
sitions for shock models with densities of 10 3 ' 5 cm -3 in units of 
erg s — 1 cm -2 arcsec -2 . The shock velocity and magnetic field b 
parameter used for each shock model are given on the top and left 
of the grid, respectively, while the model name is given in the top 
left hand corner of each box. The green (lightest) lines show the 
ISolomon etaLI 1(19871 1 12 CO J = 1 -> line strengths, the blue 
(darkest) lines show the CO shock spectra, and the red (medium) 
lines show the CO PDR spectra. The ten lowest rotational tran- 
sitions of CO are labeled in the lower right grid panel. Note how 
the shock spectra dominate over the PDR spectra for high J tran- 
sitions. 



Table 4 

Predicted ^-Integrated Intensities 



Model 




S(0) 




S(l) 




(10 


_ 1 Q _1 

±a erg s 


(10 


19 erg s 1 




cm 


arcsec ) 


cm 


2 arcsec -2 ) 


(1) 




(2) 




(3) 


n25v3bl 




1.5 




3.5 


n30v3bl 




3.0 




6.9 


n35v3bl 




5.0 
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Note. — Column 1 shows the model names 
while Columns 2 and 3 show the integrated inten- 
sities of the S(0) and S(l) lines, respectively, for 
the three shock models with significant H2 cooling. 



tional heating sources, such as cosmic-ray or photoelec- 
tric effect heating, and thus, underestimates the tem- 
peratures that would be achieved in these weak shocks. 
The combination of shock heating and cosmic-ray heat- 
ing must produce temperatures higher than that pro- 
duced by cosmic-ray heating alone. Note that in Figure 
[TJ the gas temperature is initially less than 0.1 K and 
falls below 10 K by the end of the model due to the lack 
of these extra heating terms. As such, we do not con- 
sider the models with b = 0.3 and a shock velocity of 2 
km s -1 to be valid. The lack of cosmic-ray heating does 
not significantly affect the other models because the peak 
temperatures are much larger than 10 K in these mod- 
els, indicating that shock heating is significantly stronger 



than what cosmic-ray heating would be. 

Figures [2J3] show that almost all of the emission from a 
molecular cloud in the lowest three rotational transitions 
of CO comes from unshocked gas. These figures also 
show, however, that most of the emission coming from 
the mid to high J CO lines (lines at or above J = 7 — s- 6) 
comes from shocked gas. Not only should the integrated 
intensities of these lines be higher than predicted from 
PDR models, but the excitation temperatures derived 
from the line ratios of these lines should also be higher 
than a PDR model would predict. Thus, these mid to 
high J CO lines should serve as observational diagnostics 
of turbulent energy dissipating via shocks. Since CO 
accounts for the majority of the cooling via shocks, these 
mid to high J transitions of CO should be the best tracers 
of where the majority of the energy goes in turbulent 
shocks in molecular clouds. 

Shock emission dominates at higher J transitions be- 
cause the CO in the shocked gas is warmer than the 
majority of the CO in the rest of the cloud. While the 
outer layers of a molecular cloud can be quite warm, due 
to the incident ISRF, there is little CO in these warm 
outer regions. The majority of the CO flux in the PDR 
models comes from gas that is below 20 K. Since our 
low-velocity shocks are C shocks, CO survives the shock 
and r adiates from gas at tem peratures in excess of 50 K. 

The lSolomonl^aLl (fl987l) 12 CO J = 1 -> integrated 
intensity is larger than the predicted shock-integrated in- 
tensity in all of the models, which confirms that shock 
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emission is not the major source of emission in the 1 
— > transition. While the comparison PDR models 
for the high-density shock models over p redict the 1 — » 

inte grated intensity, compared to the Solomon ct al. 
(1987)-integrated intensity, the comparison PDR models 
for the 10 25 cm -3 shock models have J = 1 — > inte- 
grated intensiti e s in re markably good agreement with the 
iSolomon et al.l (|1987j )-integrated intensity. This slight 
discrepancy be t ween some of the PDR models and the 
ISolomon et al.l (|1987l ) relation is likely due to the char- 
acteristic lin e width, depth and density of the clouds 
observed by ISolomon et al.l (|1987f ) being slightly differ- 
ent than the values used in the PDR models. 

4.2. Variation across Parameter Space 

The maximum temperature reached in the shocked 
gas increases with increasing Alfvenic Mach number and 
Mach number of the shock. All of the models with a 
shock velocity of 3 km s _1 , and thus a Mach number 
of 17, have higher maximum temperatures than all of 
the models with a shock velocity of 2 km s" 1 , corre- 
sponding to a Mach number of 12. For models with 
the same Mach number, models with higher Alfvenic 
Mach numbers, those with lower magnetic b parameters, 
have higher maximum temperatures. A larger Alfvenic 
Mach number alone, however, does not necessarily imply 
a larger maximum temperature, as the models with b 

1 and a shock velocity of 2 km s _1 have roughly the same 
maximum temperature as the models with b = 3 and a 
shock velocity of 3 km s _1 despite having almost twice 
as large Alfvenic Mach numbers. The maximum temper- 
ature reached significantly affects the CO shock profile, 
as higher temperatures excite higher rotational states, 
which leads to considerably more emission in higher J 
transitions. As described further in Section l4~4l the ef- 
fectiveness of H2 cooling is also highly sensitive to the 
maximum temperature reached in the gas. 

The fraction of energy going toward compressing the 
magnetic field is primarily dependent upon the Alfvenic 
Mach number of the shock, with more energy going into 
the magnetic field in the models with lower Alfvenic 
Mach numbers. Larger Alfvenic Mach numbers also pro- 
duce smaller cooling lengths, cooling times, and filling 
factors. 

The turbulent energy density of a molecular cloud 
is dependent upon both the density and turbulent ve- 
locity dispersion of the cloud. Thus, the models with 
higher densities and higher shock velocities have larger 
integrated intensity scaling factors (see Equation (|15p) 
and, consequently, higher integrated intensities for all 
CO transitions. 

The critical densities for all the CO line s above the J 
= 3 — > 2 line are greater than 10 5 cm -3 (iSchdier et al.1 
2005). Therefore, larger initial densities make it eas- 
ier for higher rotational states to be populated and, as 
such, more emission comes out in higher lying lines in 
the higher density models. This effect is clearly seen in 
the shock models with b = 0.1 and a shock velocity of 2 
km s , as the line with the largest integrated intensity 
shifts from the 3 — > 2 line to the 4 — > 3 and then finally to 
the 5 — > 4 line as the density increases from 10 25 cm" 3 
to 10 3 cm -3 and then to 10 3 5 cm -3 . The peak of the 
PDR spectra remains at the 3 — > 2 transition in all of 
the PDR comparison models because the same density 



is used for all of the PDR comparison models. This dif- 
ference between the shock and PDR models causes the 
transition at which shock emission becomes larger than 
PDR emission to move to slightly lower transitions as the 
density of the shock model increases. 

As described above, changes to the magnetic field 
strength, shock velocity, or initial gas temperature can 
significantly alter the shape and the scaling of the shock 
CO spectrum. None of these changes, however, alter the 
key result that shock emission dominates at mid to high 
J transitions, particularly from J = 7 — > 6 and up. 

The fraction of energy dissipated via CO cooling is 
not strongly correlated with any shock property. Only 
in the strongest shock models does the fraction of energy 
dissipated via CO weakly depend upon the initial density 
of the gas. For these strong shock models, the fraction of 
energy that is emitted via CO rotation lines increases by 
approximately 10% as the initial density increases from 
10 2 - 5 cm" 3 to 10 3 - 5 cm" 3 . 

An increase in the initial density of the gas weakly 
increases the effectiveness of gas-grain cooling, but this 
cooling term accounts for at most 4% of the shock cooling 
in any of the models. Much higher densities, densities 
closer to 10 5 cm -3 , are required before gas-grain cou- 
pling becomes reasonably efficient. Going from a density 
of 10 2 ' 5 cm" 3 to 10 3 5 cm" 3 also decreases the cooling 
lengths, cooling times, and filling factors of all of the 
models by approximately a factor of six. 

4.3. Magnetic Field Compression 

In our low velocity shock models, the energy that goes 
into compressing the magnetic field is of the same order of 
magnitude as the energy radiated away in CO rotational 
lines. It is, however, unclear where this injected mag- 
netic energy would go. A local increase in magnetic field 
strength could further drive MHD waves, which would 
subsequently shock. In this case, more shocks would be 
required in order for all of the cloud's turbulent energy 
to be dissipated by CO cooling and our predicted line 
integrated intensities would have to be increased by a 
factor of two. 

Alternatively, this magnetic energy may slowly leak 
out of the cloud via magn etic coupling with the external 
medium, as de scribed by Elmegreenl (|1985l ) and seen in 
simulations by lEnd (|2002f ). This magnetic energy may 
also be dissipated on small scales via a process such as 
ambipolar diffusion. 

4.4. H 2 Lines 

Molecular hydrogen does not have a permanent dipole 
moment and thus, radiates through weak quadrupole 
transitions (AJ = 2). Furthermore, since hydrogen is 
so light, the rotational energy levels of molecular hydro- 
gen are relatively widely spaced. These two effects make 
H2 rotational emission highly temperature sensitive and 
temperatures in excess of 100 K are required for signifi- 
cant emission. This temperature sensitivity can be seen 
in the shock models, as H2 emission is essentially negligi- 
ble in all but the strongest shock models, which are the 
only models where the maximum temperature exceeds 
100 K. In these models, the S(l) and S(0) H 2 lines have 
comparable integrated intensities to the CO J = 7 — >• 6 
line, although the H2 lines are relatively stronger in the 
lower density models. 
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The lack of H2 rotational emission from cool gas means 
that, aside from shocked gas, the only significant source 
of H2 emission in a molecular cloud is the thin outer edge 
of the cloud's PDR where the temperature is high and 
H 2 is not rapidly photodissociated. Thus, H 2 rotational 
emission could be a very useful tracer of shocked gas in a 
molecular cloud. In particular, while the predicted S(0) 
shock-integrated intensities are on the order of the S(0) 
PDR-intcgratcd intensity, the S(l) shock- integrated in- 
tensities range from being five times to over fifteen times 
larger than the S(l) PDR- integrated intensity. The ra- 
tio of S(1)/S(0) at the 150 K maximum temperature 
of the strongest shocks, approximately 2.3, is also sig- 
nificantly different fro m the ratio of these lines in the 
iKaufman etldl ([2006) comparison PDR models, roughly 
0.3. 

H2 cooling may be even more significant in gas that is 
lacking in gas phase CO, as this shocked gas is likely to 
reach higher temperatures with the effectiveness of CO 
cooling reduced. This increase in H2 shock emission in 
CO sparse gas, such as th e "dark gas" in the periphery 
of a molecular cloud (e.g.. IWolfire et aTlEJOK):. may nat- 
urally produce a limb brightening effect for H2 rotational 
emission in molecular clouds, as observed in Taurus by 
IGoldsmith et~aTl (l20Toh . 

In all of our shock models, cooling from vibrational 
transitions of H2 is negligible because temperatures on 
the order of a few thousand Kelvin are required to excite 
higher energy vibrational states (jKaufman fc Neufeldl 
Il996al) . 

4.5. Other Shock Tracers 

Water is kno wn to be an effective coola nt in high- 
velocity shocks (jKaufman fc Neufeldl Il996al ) but water 
cooling is negligible in all of our shock models. This 
is because water is only forme d in gas with a tempera- 
ture of a few hundred Kelvin (|Elitzur fc de Jong] 119781 : 
lElitzur fc WatsonH l978) and is only efficiently liberated 
from du st grains, due to sp uttering, in 15 km s _1 or faster 
shocks (|Draine et al.lll983l) . Low- velocity shocks are also 
very ineffective at heating dust grains, meaning that 
thermal sublimation of water off of dust g rains is com- 
plete ly negligible in our low- velocity shocks (jDraine et al.1 
fl983h . 

At the low densities of our shock models, the interac- 
tion timescale of gas with dust grains is long compared 
to the cooling time such that only a few percent of the 
energy is ever liberated via gas-grain interactions in the 
models. Densities closer to 10 5 cm -3 are needed before 
gas-grain coupling becomes effective. 

While CO line radiation and the compression of mag- 
netic fields are the dominant coolants in our low-velocity 
shocks, other molecular lines, which have not been in- 
cluded in our shock models, may still be valuable tracers 
of shocked gas. In particular, molecular lines that are 
sensitive to increased temperature could provide shock 
tracers in different wavelength regimes. Molecular tran- 
sitions that are sensitive to density may also be useful 
shock tracers, but are less likely to be as useful as temper- 
ature sensitive lines since the maximum densities reached 
in the shocked gas are less than the gas density of prestel- 
lar cores. 

4.6. Additional Caveats 



In scaling the shock models to predict the total 
strength of the shock emission from an entire molecular 
cloud, it was assumed that all of the turbulent energy of 
the cloud is dissipated at one particular shock strength. 
In reality, energy will be dissipated through a variety of 
different strength shocks, either due to different shock ve- 
locities or different strengths of the magnetic field paral- 
lcl to the shock front. Furthermore. iSmith et all ([2000a) 
show that high- velocity shock s dominate energy d issipa- 
tion in driven turbulence while ISmith et all ()2000bf ) show 
that if the turbulence is decaying, low-velocity shocks 
dissipate most of the energy. If energy is dissipated 
through lower velocity shocks, then our calculated line- 
integrated intensities will overestimate the actual emis- 
sion from molecular clouds in lines with higher excita- 
tion temperatures. If turbulence is driven at scales much 
smaller than the size of the cloud, however, then our line 
strengths should be increased by a factor of k _1 (see Sec- 
tion [221 for a discussion on how k relates to the driving 
scale of turbulence). 

Another factor of two uncertainty in the shock- 
integrated intensities comes from the assumption that 
the energy in magnetic field fluctuations in the cloud is 
negligible. While some MHD simulations have shown 
that the turbulent kinetic energy dominates over the 
energy in magnetic field flu ctuations, particularly for 
smaller initial magn e tic fields dPadoan fc Nordlundll 19991 : 
iPadoan et afl [2000t iHeitsch et all 120011 ). other simula- 
tions find that these magnetic field fluctuations have en- 
ergy on the order of the kinetic energy of the turbulence 
(jStone et al. 1998: IQstriker et al.lf200lh . If the energy in 
magnetic field fluctuations is on a par with the kinetic en- 
ergy of turbulence, the line-integrated intensities would 
have to be scaled up by a factor of two to account for the 
dissipation of this additional energy. 

In scaling up the shock models to estimate the total in- 
tegrated intensities from an entire cloud, the integrated 
intensity of every line in a particular model has been 
multiplied by the same factor and no optical depth ef- 
fects have been taken into account. These optical depth 
effects, however, should not be significant for the higher 
rotational transitions of CO, where the CO transitions 
are effective at tracing shock emission, as the PDR mod- 
els indicate that the CO lines are only optically thick 
up to, and including, the J = 5 — > 4 transition. As for 
shock emission in the lower rotational transitions of CO, 
this emission is likely to be absorbed by the ambient gas 
and thereby will serve as a heating source for the ambi- 
ent gas. Note, however, that the expected CO-integrated 
intensity at these low transitions is much less than the 
PDR-integrated intensity, which indicates that this extra 
shock heating will have only a very minor effect on the 
temperature, and thus the spectrum, of the ambient gas. 

In deriving the total turbulent energy dissipation rate 
of a molecular cloud, it was assu med that turbulenc e de- 
cays in a cross i ng tim e. Recently, iBasu et ail (|2009fl and 
iBasu fc Dappl (|2010l) discovered a long-lived magnetic- 
tension-driven mode in their thin disk simulations of 
flattened molecular clouds, arising from interactions be- 
tween the disk and an external magnetic field, which was 
able to preserve a significant fraction of the turbulent 
energy of the cloud for much longer than the crossing 
time. The existence of such a long-lived MHD mode 
would significantly reduce the required energy dissipation 
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rate of a molecular cloud and, therefore, our predicted 
line-integrated intensities as well. This long-lived mode, 
however, has not been noticed in any further simulations, 
including the t hree-dimensiona l simu lations of collapsing 
cores done by Kudoh fe Basu | (1201 If), who d id loo k for 
this particular MHD mode. iKudoh fe Basul (|2011[ ) sug- 
gest that the absence of this mode in their simulations 
may be due to the very small density contrast between 
the disks and surrounding gas in their simulations. 

The line ratios from these shock models are indepen- 
dent of the scaling of the lines and thus, are not affected 
by the above uncertainties. The line ratios of the shock 
models are also independent of the assumption of spher- 
ical geometry. 

The turbulent energy of molecular clouds may also not 
be dissipate d completely through shocks. iSmith et al.1 
(|2000b[ ) and lStone et al.l (jl998| ) find, in their simulations 
of turbulence, that only 50% of the turbulent energy is 
dissipated through artificial viscosity, due to the presence 
of shock fronts, while the other 50% is dissipated through 
numerical viscosity, representative of small-scale dissipa- 
tion distributed relatively uniformly across the cloud. It 
is possible that some of this uniformly dissipated energy 
should have been dissipated in weak shocks or vortices 
that were not resolved in these simulations and thus, we 
consider that 50% is only a lower limit for the fraction of 
energy dissipated in intermittent structures (i.e., shocks). 

An alternative model for turbulent dissipation, where 
energy is dissipated thr ough magnetized vo rtices, has 
also been pu t forw ard (jGodard et al.l [2009). In the 
iGodard et al.1 (|2009f ) turbulent dissipation region (TDR) 
model, small vortices on the order of a few tens of AU 
heat gas to temperatures of nearly 1000 K via ion- neutral 
friction. The spectral signature of such a TDR model 
should be significantly different from our low-velocity 
shock models because the TDR model produces temper- 
atures much higher than what can be obtained from slow 
shocks. 

All of the shock models have been run with standard, 
roughly solar, metal abundances. In lower metallicity 
clouds, the abundance of CO will be reduced. The gas 
phase CO abundance will also be lower at higher den- 
sities, due to fre eze out of CO onto dust grains (e.g., 
iGoldsmithl |2001[ ). and in the outer layers of molecu- 
lar clouds, where CO is readily photodissociated (e.g., 
iWolfire et al.l l2~010"). The reduction of gas phase CO may 
lead to higher post shock temperatures, which would 
change the resulting CO spectrum and could affect which 
cooling mechanism is dominant. H2 line cooling and the 
deposition of energy into magnetic fields may also be 
more important for dissipating kinetic energy in CO poor 
gas. Further low-metallicity shock models, however, are 
needed to confirm this. 

The CO spectrum predicted from the PDR model is 
dependent upon the chemistry put into the models. In 
particular, any chemical effects which would increase the 
temperature in the outer CO layers, such as an increase 
in polycyclic aromatic hydrocarbon heating, would shift 
the PDR spectrum toward higher rotational transitions. 
The PDR models used also have a relatively low ISRF 
of 3 Go- In active, high-mass star- forming regions, the 
ambient ISRF is likely to be much higher, due to the 
presence of previously formed, massive, young stars and 
thus, the PDR emission from these regions is likely to be 



much stronger and peaked toward much higher J tran- 
sitions. As such, our model comparisons should only be 
used for relatively quiescent, low-mass star-forming re- 
gions in which the ambient ISRF is relatively low, rather 
than in strong UV environments such as the Orion Bar. 

4.7. Observational Potential 

The total energy dissipated by shocks cannot be di- 
rectly observed because some of the turbulent energy is 
not radiated away, but rather goes into increasing mag- 
netic field strengths. Furthermore, many of the low lying 
CO lines are not readily observable, because these tran- 
sitions are dominated by emission from unshocked gas. 
This ambient gas emission, however, does drop rapidly 
at higher J numbers and the CO J = 6 -> 5 (691.47308 
GHz) and 7 — >• 6 (806.651806 GHz) transitions are dom- 
inated by shock emission in most of our models. Thus, 
these higher rotational CO transitions act as shock trac- 
ers and by fitting shock models to the observed strengths 
of multiple high J CO transitions, the total shock lumi- 
nosity of a cloud can be estimated. 

The Herschel Space Observatory's Heterodyne Instru- 
ment for the Far Infrared (HIFI) and Spectral and Pho- 
tometric Imaging Receiver (SPIRE) both have the neces- 
sary wavelength coverage and sensitivity to be able to de- 
tect the CO J = 5 — ^ 4, 6 — s> 5, and 7 — » 6 lines, if they are 
as bright as we predict in our shock models. These two 
Herschel instruments also cover the J = 8 — > 7 wavelength 
but our predicted line strengths for this transition are at 
the limits of what could be detected within a few hours 
of observing time. From the ground, the James Clerk 
Maxwell Telescope's (JCMT) receiver W is capable of 
observing at the wavelength of the CO J = 6 — > 5 transi- 
tion while the Atacama Pathfinder Experiment (APEX) 
and the Caltech Submillimctcr Observatory (CSO) have 
instruments that operate in the appropriate wavelength 
regimes to detect both the 6 — > 5 and 7 — » 6 lines. The 
CO J = 6 — s- 5 and 7—^6 lines also lie within bands 
9 and 10, respectively, of the Atacama Large Millimeter 
Array (ALMA). ALMA, with its superb resolution, may 
resolve individual shock fronts and thus, provide infor- 
mation regarding the properties of individual shocks. 

Some care must be taken when choosing a location in 
a molecular cloud to observe, as there are other sources 
of high J CO line emission that have not included in 
the comparison PDR models. High-velocity protostel- 
lar outflows will generate large shocks that can eas- 
ily produce temperatures of hundreds of Kelvin. Such 
strongly shoc ked gas will radiate signifi cantly in high 
J lines (e.s... Kaufman fe Neufeldlll996blial ). Embedded 
high-mass stars will also significantly heat nearby gas and 
lead to high J line emission. Thus, while turbulent dissi- 
pation should occur in active, high-mass star-forming re- 
gions, the spectral signatures of low- velocity, turbulence- 
induced shocks may only be readily detectable in qui- 
escent regions of low-mass star-forming molecular com- 
plexes. Water emission can be used as a discrimi- 
nant between low-velocity, turbulent shocks and the 
stronger shocks produced by outflows, as little water 
emission is predicted from our models while stronger 
shock s are expected to cool sign ificantly via water lines 
fe.g.. Kaufman fc Neufeldlll996aj) . 

The temperature sensitivity of the S(0) and S(l) lines 
of H 2 , at 28.2 /zm and 17.0 /xm, respectively, makes these 
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lines potential shock tracers. The atmospheric transmis- 
sion at 28.2 /Ltm is, however, very poor, and thus, it is 
extremely difficult to observe the S(0) line with ground- 
based facilities. The expected weakness of these two lines 
makes observations of H 2 emission even more problem- 
atic. 

5. GLOBAL HEATING 

5.1. Cosmic Rays 

In the cold, well-shielded, central regions of molec- 
ular clouds, cosmic-ray heat ing is believed t o be the 
dominant heating term (e.g.. iGoldsmithl [20011 ). but the 
general cosmic-ray ionization rate in the galaxy is not 
particularly well known. Estimates for the cosmic- 
ray ionization rate i n dense gas vary from 10 ~ 16 s _1 
(|Caselli et all [l99l ILTsitl [20031 Ibotv et all l200l to 
10~ 18 s" 1 (ICaselli et all [TIM [20021: [Flower et al.l I2007t 
IHezareh et al.l I2008T) but commonly lie between 1 x 
IP" 17 s" 1 and 5 x 1Q- 17 s " 1 (IWilliams et all ~ 
van der Tak fc van Dishoeckl l200u iDotv et aTT 



1998 



2002 



Wakelam et alJl2005HBergin et al . 2006; Maret fc Berginl 
20071: IGoicoechea et all I2009D . While the spread in 



cosmic-ray ionization rat es may simply be due to spa- 
tially varying rates (e.g.. Ivan der Tak et al.1 12006). some 
of the scatter is likely due to uncertainties in the chem- 
ical models used to determine the ionizati on rates, 
as ch anges in known reaction rates (e.g., IDalgarnq 
I2006D and the inclusion of non-equilibrium chemistry 
( Lintott k, Rawlings 2006) have changed the results of 
ionization rate calculations. 

For this paper, a cosmic-ray ionization rate range of 
1 x 10 -17 s _1 to 5 x 10 -17 s" 1 is adopted and, since each 
cosmic-r ay ionization depo sits approximately 20 eV into 
the gas ([Goldsmith! I2001f) . the volume heating rate by 
cosmic-rays is taken to be 

r cr = 0.3 to 1.6 x 10~ 27 f — erg s~ x cm" 3 . (19) 

\cm~ A J 

This cosmic-ray heating rate is shown in Figure [5] 

5.2. Turbulent Heating 

If the turbulent energy of a cloud is not dissipated 
through localized shocks, but rather is dissipated rela- 
tively uniformly across the cloud, then this turbulent en- 
ergy dissipation will act more as a general heating mech- 
anism for the cloud and will yield the heating rate given 
by Equation For the sole purpose of comparing the 
turbulent energy dissipation rate to the cosmic-ray heat- 
ing rate, b oth the size- velocity a nd density-size relation- 
ships from lSolomon et al.l ([1987f ) are used to rewrite the 
turbulent energy dissipation rate in terms of only n and 
the gas density. Using these two scaling relations, the 
turbulent energy heating rate can be written as 



turb 



= 5.86 x 10 



-25 -1 



10 3 ( 



0.5 



erg cm 



-3 



( 2 °) 

where e is the fraction of the dissipated turbulent energy 
th at acts as a g l obal he ating mechanism. The s imulations 
of ISmith et all (|2000bD and lStone et~aTI (|1998l ) both sug- 
gest that roughly half of the turbulent energy of a cloud 
may be dissipated uniformly across a cloud instead of in 
localized shocks and thus, suggest that e = 0.5 (see Sec- 
tion [4]6] for further discussion on e). Furthermore, the 



magnetic energy injected by shocks may also decay and 
lead to a global heating of the cloud. Figure [5] shows 
the turbulent heating rates, as a function of density, for 
e = 1 and e = 0.5. As before, a fixed k value of 1 is used. 

5.3. Ambipolar Diffusion 

One potentially important energy dissipation mecha- 
nism that is usually not include d in simulations (e.g., 
iStone et a,l.lll998t iMac Lowl[l999t) i s ambipolar di f fusion 
From their numerical simulations, [Padoan et al.l ([2000D 
find that the ambipolar heating rate in well shielded 
molecular gas is given by 

r - \( Ma \ 2 ( 



( <\B\> 
V 10 \iG 

io- 24 



< n > 
,320 cm- 3 



erg cm 3 s , 



(21) 



where < \B\ > is the volume-averaged magnetic field 
strength, Ma is the Alfvenic Mach number of the tur- 
bulence, and < n > is the volume-averaged number den- 
sity. Applying Larson's laws and a density-magnetic field 
strength scaling relation with k = 0.5 converts this equa- 
tion to the form 



ambi 



= 2.6 x 10" 



10 3 cm 



-0.5 



erg cm 3 s . 



( 22 ) 

Figure [5] shows this ambipolar diffusion heating rate, as 
a function of density, evaluated for the two b parameters 
previously used for our shock models, b = 0.1 and 0.3. 
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Figure 5. Various heating rates for well-shielded molecular gas. 
The shaded (green) region shows the range of cosmic-ray heating 
rates, the dark (purple) solid line shows the total shock turbulent 
energy dissipation rate, the light (blue) solid line shows 50% of 
the shock turbulent dissipation rate, and the dark (red) and light 
(yellow) dotted lines show the ambipolar diffusion heating rates for 
b values of 0.3 and 0.1, respectively. Please see the online version 
for a color version of this figure. 



5.4. Discussion 

For gas with a density larger than 10 3 cm~ 3 , the 
cosmic-ray heating rate is the dominant heating term. 
The turbulent energy dissipation rate is comparable to 
the cosmic-ray heating rate around a density of 10 3 cm~ 3 
and becomes larger than the cosmic-ray heating rate at 
lo wer densiti es. 

IGoldsmithl (pOOlh examines the thermal balance of 
molecular clouds and finds that the ambient tempera- 
ture of well-shielded gas, at a density of 10 3 cm -3 , is 
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10 K. We took the iGoldsmithl (|200lD prescriptions for 
cooling and heating and increased the cosmic-ray heat- 
ing by a factor of two to reproduce the effect of having 
a turbulent energy dissipation heating rate equivalent to 
the cosmic-ray heating rate. With this increased heating 
rate, we find a slightly high er equilibrium temperature of 
13 K. Pan & Padoan ( 2009) present a more detailed ther- 
mal balance model that also includes a turbulent heating 
term and, for a 1 pc sized cloud, find similar gas tempera- 
tures of 13-17 K, depending upon the cosmic-ray heating 
rate used. 

This small change in temperature is well within the in- 
trinsic scatte r of observed gas tempe ratures in molecular 
clouds (e.g.. iBergin fc Tafallal [20071. Since the cosmic- 
ray heating rate is also uncertain by at least a factor 
of two, the finding of previous thermal balance stud- 
ies that the ambient temperatures of molecular clouds 
are roughly con sistent with heating by cosmic-ray ioniza- 
tion alone (e.g.. IBergin et aLll2006D is not in conflict with 
the presence of heating from turbulent dissipation at the 
rate calculated above. The above agreement between ob- 
servations and models does, however, constrain the tur- 
bulent heating rate to not be significantly greater than 
estimated above. This implies that k c annot be much 
less th a n one, similar to the findings of iBasu fc Murall 
(|2001f) . iPadoan et all (|2000l) also note that if turbulent 
heating is significant in a molecular cloud, then a posi- 
tive temperatu re- velocity d i spersi on relation is expected, 
as observed bv lJiiina et al.l ()1999h . 

The significance of ambipolar diffusion to the thermal 
balance of a cloud is highly dependent upon the strength 
of the magnetic field. For a strong field, b = 0.3, am- 
bipolar diffusion should be the dominant heating pro- 
cess in well-shielded gas with an average density of 100 
cm~ 3 , and should be comparable to other heating pro- 
cesses at a density of 10 3 cm -3 . Gas at densities of 100 
cm~ 3 , however, may not necessarily be well shielded and 
may instead be dominated by the ISRF. If magnetic field 
strengths are slightly lower, corresponding to b = 0.1, 
the ambipolar diffusion rate is negligible for all densities 
greater than or equal to 100 cm -3 . It should also be 
noted that the ambipolar diffusion heating rate becomes 
larger at lower densities, unlike the turbulent and cosmic- 
ray heating rates, because the typical collision speed be- 
tween ions and neutrals increases with decreasing density. 

Shocks will cause both the density and magnetic field 
strength to increase. The change in the density will be 
proportional to the change in the magnetic field strength, 
given the flux freezing approximation made in the shock 
models. The Alfven speed in the gas will thus also in- 
crease by a factor proportional to the square root of the 
change in the magnetic field strength. If the magnetic 
field strength is locally increased by a factor of q by a 
shock and the velocity dispersion of a cloud remains rel- 
atively unchanged, then the ambipolar diffusion heating 
rate in the shocked gas will be 

r ambi. shocked — Q Pambi.Oi (^^) 

where r a mfn.o is the ambipolar diffusion dissipation rate 
in the unshocked gas. Thus, the postshock gas will have 
its ambipolar diffusion heating rate increased. 

Table [3] gives the factors by which the magnetic field 
strength and the ambipolar diffusion heating rate in- 



crease in each of the shock models. The magnetic field 
increases in strength by a factor of 4-22 and the ambipo- 
lar diffusion heating rate increases by one to two orders 
of magnitude. Even for the weakest magnetic field mod- 
els, the ambipolar diffusion heating rate in the shocked 
gas is greater than the rate at which energy is injected 
into the magnetic field through the shock and thus, this 
enhanced ambipolar diffusion rate may provide an ef- 
ficient mechanism for dissipating the injected magnetic 
energy. Furthermore, this enhanced ambipolar diffusion 
rate should add an extra heating source to the shocked 
gas, thereby increasing the CO flux that will emerge at 
higher J transitions. 

The above discussion regarding heating rates is only 
relevant for the well-shielded centers of molecular clouds, 
as photoelectric heating due to the ISRF will dominate 
the h eating in the outer P DR zones of molecular clouds 
(e.g.. iKaufman et al.lll999T ). The heating rates shown in 
Figure[S]were also de rive d using Larson's laws, which, as 
discussed in Section 12.11 have recently been called into 
question. 

6. CONCLUSIONS 

We have run models of MH D, C-type shocks, based 
on IKaufman fc Neu fcld ( 1 996bh . for shock velocities of 2 
and 3 km s _1 and initial densities between 10 2 5 and 10 3 5 
cm~ 3 . CO is found to be the dominant molecular coolant 
with 40%-80% of the shock energy being emitted in CO 
rotational lines. All other line cooling processes are neg- 
ligible, except for H2 line cooling in the models with the 
very strongest shocks, in which H2 cooling accounts for 
5%-20% of the total shock cooling. Between 20% and 
60% of the shock energy also goes into compressing the 
magnetic field. 

The expected CO spectrum from each of the shock 
models has bee n calcu lated and PDR models, based on 
IKaufman et al.l (|1999|) . were used to determine the ex- 
pected contribution of CO emission from not only the 
cold, well-shielded interior of a molecular cloud, but also 
from the warm outer layers of the cloud. The PDR emis- 
sion dominates for low J transitions of CO but the shock 
emission is larger at mid to high J transitions. In all 
models the shock emission is larger than the PDR emis- 
sion in the J = 7 — > 6 transition and the shock emission 
can dominate as low as the J = 5 — > 4 transition, de- 
pending upon the shock model. The J = 6 — > 5 and 7 — > 
6 should serve as shock tracers and should be detectable 
with current observational facilities. 

The turbulent energy dissipation rate is larger than 
the cosmic-ray heating rate for densities less than 10 3 
cm~ 3 . The presence of such an additional heating term 
is, however, still consistent with previous thermal balance 
studies, given the uncertainty in the cosmic-ray heating 
rate and the range of observed molecular cloud tempera- 
tures. The ambipolar diffusion heating rate is negligible 
at high densities and low-magnetic field strengths but 
can be dominant at lower densities and higher magnetic 
field strengths. Ambipolar diffusion is also enhanced in 
the shocked gas and may provide a mechanism for the 
dissipation of energy injected into the magnetic field by 
a shock. 
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role of magnetic fields in molecular clouds as well as our 
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